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Abstract 

Superfluid hydrodynamics affects the spin-evolution of mature neutron stars, and may be key to 
explaining timing irregularities such as pulsar glitches. However, most models for this phenomenon 
exclude the global instability required to trigger the event. In this paper we discuss a mechanism 
that may fill this gap. We establish that small scale inertial r-modes become unstable in a superfluid 
neutron star that exhibits a rotational lag, expected to build up due to vortex pinning as the 
star spins down. Somewhat counterintuitively, this instability arises due to the (under normal 
circumstances dissipative) vortex-mediated mutual friction. We explore the nature of the superfluid 
instability for a simple incompressible model, allowing for entrainment coupling between the two 
fluid components. Our results recover a previously discussed dynamical instability in systems where 
the two components are strongly coupled. In addition, we demonstrate for the first time that the 
system is secularly unstable (with a growth time that scales with the mutual friction) throughout 
much of parameter space. Interestingly, large scale r-modes are also affected by this new aspect 
of the instability. We analyse the damping effect of shear viscosity, which should be particularly 
efficient at small scales, arguing that it will not be sufficient to completely suppress the instability 
in astrophysical systems. 



I. INTRODUCTION 

Since the first observation of a glitch in the spindown of the Vela Pulsar in 1969 jl, 2], 
there have been a large number of similar events observed in other pulsars {3, 4]. Several 
competing, in some cases complimentary, mechanisms have been suggested as explanation 
for these occurrences (5 9 1 - The two most widely classes models relate either to changes in 
the star's elastic crust or the dynamics of the superfluid neutrons that are present both in 
the core and the inner crust. The first set of models involve fractures in the crust, leading 



to changes in the moment of inertia of the star [lOj. The second set involves the unpinning 
of vortices associated with a superfluid component from the star's inner crust 11] . The 
unpinning event is attributed to the build-up of "differential rotation" between the two 
components (elastic crust and interpenetrating superfluid). It is the second of these two 
mechanisms with which we concern ourselves in this paper. 

We consider a promising mechanism for triggering the observed events; a superfluid in- 
stability present in systems with "differential rotation" , like the lag that plays a central ro 



in all vortex-based glitch models. The basic mechanism has already been discussed in [12], 
where it was argued that the instability may trigger large-scale vortex unpinning leading to 
the observed events. The first part of the argument was a demonstration that unstable r- 
modes could be generated by the difference in rotation of the charged particles (protons and 
electrons) and the superfluid neutrons in the star. In particular, it was shown in 12] that 
such unstable modes exist for a strongly coupled system. The second part of the argument 
consisted of a demonstration that viscosity would not suppress this instability completely, 
but would allow growing modes to exist for a range of wavelengths once a critical rotational 
lag was reached. The third part of the argument placed the mechanism in an astrophysical 
context by comparing the predictions of the model to observational data. 

In this paper we consider this new r-mode instability in more detail. We extend the 
analysis of [12] beyond the strong-drag limit, and show that dynamically unstable modes 
(with a growth time similar to the star's rotation rate) exist for a wide range of parameter 
values. In addition, we discover that the r-modes also suffer a secular instability (with 
growth time proportional to the mutual friction) throughout much of parameter space. This 
aspect is new. It is particularly interesting as it affects also the global scale r-modes, while 
the dynamical instability is restricted to small scale modes. As in the previous work, we 
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imit ourselves to discussion of the non-relativistic case. Again following the approach of 



12] , we introduce shear viscosity and demonstrate that the instability is not completely 



suppressed by the inclusion of the associated damping. 

The main purpose of this investigation is to establish the robustness of the superffuid 
instability, highlight the key ingredients that lead to its presence and set the stage for 
more detailed studies of this mechanism. The instability that we consider may be generic 
(belonging to the class of two-stream instabilities that are known to operate in a wide variety 
of physical systems [l3|), but at this point it is not clear how it is affected by other aspects 
of neutron star physics. In particular, future work needs to extend the analysis to account 
for the crust elasticity and the presence of the star's magnetic field. 



II. SETTING THE STAGE 



We consider the conditions that are expected to prevail in the outer core of a mature neu- 
tron star, where a neutron superfluid is thought to coexist with a proton superconductor. 
A key qualitative aspect of this system is that it allows the two "fluids" to flow "indepen- 
dently" , leading to dynamics that is well represented by a two-fluid model. In order to keep 
the analysis manageable, we assume that both fluid components are charge neutral. This 
allows us to neglect (in this first proof-of-principle analysis) the star's magnetic field. This 
is obviously a severe simplification, but to analyse the full problem would be difficult at this 
point. The key aspect that we add to previous studies in this problem area is the relative 
rotation between the two components, expected to build up as the observed component 
spins down due braking associated with the exterior magnetic field. We are interested in 
the global oscillations of a star that exhibits this kind of "differential" rotation. Despite 
this being a key aspect of the generally accepted "explanation" for observed pulsar glitches, 
there have not yet been any detailed studies of the dynamics of such a system. In this sense, 
the present work provides an important step towards realism. 

The main question that motivates us stems from the discovery that two-stream instabil- 

n 

ities are generic in the two-fluid model [14|. Given this, it is relevant to ask whether the 
rotational lag that builds up in a mature neutron star may lead to such instabilities and, if 
so, what the observational consequences may be. Conversely, given that we do not yet know 
what the mechanism that triggers the observed glitches may be, but that it ought to be 
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some kind of "global" instability, it is interesting to ask whether the two-stream instability 
may be relevant in this context. 

We consider the linear perturbations of a system that exhibits a rotational lag. Because 
of this set-up, it is natural to work in the inertial frame (rather than choosing one of the 
rotating frames). It is also natural, given the complex nature of the perturbed velocity 
fields etcetera to carry out the analysis in a coordinate basis. This means that vectors, 
like the velocity, are expressed in terms of their components, v l (say), and a distinction is 
made between co- and contravariant objects, with the former following from the latter as 
Vi — 9ij v ^ where is the flat three-dimensional metric. This description of the problem is 
also advantageous since it involves the use of the covariant derivative V, associated with the 
given metric, which automatically encodes the scale factors associated with the curvilinear 
coordinates (r, 9, ip) that are appropriate for the problem. 

We consider a system with two interpenetrating fluids, labeled n and p (from now on), 
which are assumed to rotate rigidly in such a way that 

4 = 04 , x = n,p (1) 

where e 1 represents the azimuthal basis vector (and e l 9 later the polar one). The angular 
velocities, Q x , are assumed to be constant, but we allow for differential rotation, i.e. Q l n 7^ Q l . 
To keep the analysis tractable, we assume that the two rotation axes coincide, and use 

fij, = Q i , and Qi = (1 + A)fi* (2) 

Making contact with pulsar observations, Q l would be the observed (angular) rotation fre- 
quency (i.e. that of the crust) while Q l n represents the rotation of the unseen (interior) 
superfluid component. As the superfluid is expected to lag behind as the crust spins down, 
one would typically expect A to be small (of the order of 10 -4 at the time of a glitch Q) 
and positive. 

Key to understanding the dynamics of glitching pulsars is the appreciation that the 
neutrons may form Cooper pairs, and hence act as a superfluid. The upshot of this is that 
bulk rotation can only be achieved by the formation of an array of quantized vortices. As 
discussed in, for example, [ijj] the quantization condition is imposed on the momentum that 
is conjugate to the velocity v l n . This momentum is given by 

P- = m p « + e n wD (3) 
4 



where m p is the nucleoli mass (we ignore the small difference between the bare neutron and 
proton masses) , e n represents the entrainment between neutrons and protons (more of which 
later) and the relative velocity is 

<4x = 4 - < > y + x ( 4 ) 



Given this, and the relevant quantization condition (see 15] for discussion), the local vortex 
density (per unit area), n v , follows from 



n 



,k* = e^ fc V> fc n + = + £n(^ p - K)] = 2[1 + (1 - e n )A]n* (5) 



where k 1 is the vector aligned with the vortex array with magnitude k = h/2m p (h is Planck's 
constant). Here we have introduced yet another simplifying assumption; we have taken the 
fluids to be incompressible which implies that the entrainment coefficient may (at least in 
the small A case that we are focussing on) be taken to be constant. It is also worth noting 
that p n e Q = p p E p where p x = m p n x are the respective mass densities. 

Let us now consider the linear perturbations (represented by 5) of this kind of configu- 
ration. Assuming that the flow is incompressible also at this level (incidentally not a bad 
approximation for the inertial flows that we will consider) we have, first of all, 

= (6) 

Meanwhile, the perturbed momenta 

5rf = 5vf + eJwT (7) 
are governed by the Euler equations; 

£f = tu5p* + 8viV jP * + viVjSp* + e^SwfVivi + wfVM) + V^ x = 6ft (8) 

where 

^ x = $ + j2 x (9) 

combines the gravitational potential $ and the chemical potential jl x = /i x / m p for each fluid. 
We have assumed that the time-dependence is harmonic, cx exp(iut), as we are interested 
in the oscillation modes of the system. Later, we will make extensive use of the equations 
that govern the perturbed vorticity. Specifically, we will work with the quantity 

W x = e^V^t (10) 



The right-hand side of jSJ) can be used to account for any "external" forces acting on 
the fluids. It can also be used to incorporate interactions between them. The main such 
interaction force is the mutual friction that arises due to the presence of the quantized 
vortices [if], [l^] ■ The unperturbed force takes the form [3] 

ft = ^B'n v e ljk ^w k xy + ^Bn v e ijk kh klm mw^ (11) 

Px Px 

where E and £>' are coefficients to be determined from microphysics. In the standard sce- 
nario these coefficients can be expressed in terms of a single resistivity 7Z, associated with 



scattering off of the vortices, such that 



18j 



and 

The weak- and strong coupling limits discussed later correspond to, respectively, 72. — ?• and 
1Z — > oo. It is worth noting that the system can build up a differential rotation lag only as 
long as some additional force (like vortex pinning) prevents the mutual friction from acting. 
The presence of such a force is implicitly assumed in the following discussion. 
Perturbing ( ITT]) , we get 



5f* = —B' [8(n v K j )e ijk w k y + n v K j e ijk 8w^\ 

Px 

P n K?, Mm 



To evaluate this we need 



and 



+ ^Be ljk e klm [SfaKfiffw* + nvKiivftSk? + K j 5w%)] (14) 

Px 



<5(n^) = e ljk VjSp n k (15) 



5k l = — (5} - k%) e jlm Vi5p^ (16) 



where k l is a unit vector aligned with k\ together with, cf. (J5J, 

n v K = 2[1 + (1 - e n )A]fi (17) 

In the perturbation equations discussed in the next section, we will not work with the 
momentum equations (jHJ directly. Rather, we use two combinations that represent the total 
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(perturbed) momentum and the difference. It is well established that these combinations 
isolate the two dynamical degrees of freedomin an "uncoup.ed" two-fluid system EL Hence, 
this decomposition is often used in studies of oscillating superfluid neutron stars [20| . This 
means that we need; 

P*8f? + Pv fif? = (18) 

and 



8ft -8ft =--5 ft 



(19) 



where we have introduced the proton fraction x p = p p /(p n + Pp)- Only the second degree of 
freedom is explicitly affected by the mutual friction. 

Finally, we want to account for the presence of shear viscosity (in the superfluid case 
mainly due to electron-electron scattering). In the incompressible case, this means that we 
add a force to (the right-hand side of) the proton equation of form; 



(20) 



where v ee is the kinematic viscosity coefficient. We ignore bulk viscosity for two reasons. 
First of all, glitching pulsars are cold enough that shear viscosity should be the dominant 



damping mechanism. Secondly, the particular c 



efficiently damped by bulk viscosity, anyway [32 1 . 



ass of fluid motion that we consider is not 



III. THE AXIAL PERTURBATION EQUATIONS 

Despite the various simplifying assumptions, the general perturbation problem is challeng- 
ing. It is well-known that, even in the slow-rotation limit where the star remains spherical, 



there exists a class of inertial modes [22| which require the coupling of many (spherical har- 
monics) multipoles for their description. We are, however, not going to attempt to solve the 
general problem. Instead we will ask a very specific question; How does the presence of the 
rotational lag affect the r-modes of the system? This question is relevant for a number of 
reasons, perhaps the most important being related to the fact that the large scale r-modes 
may be driven unstable by the emission of gravitational radiation [32j. From a practical 
point-of-view, it is natural to focus on the r-modes since they are associated with partic- 
ularly simple velocity fields. The hope would be that the corresponding problem remains 
tractable even when we add the rotational lag. 



The r-modes are a subclass of inertial modes (restored by the Coriolis force) that are 
purely axial/toroidal to leading order. This means that the perturbed velocities take the 
form 
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r 2 sin 6 1 r 2 sin 6> ^ 

where YJ m (0, <p) are the standard spherical harmonics, U^(r) are the mode amplitudes, and 
the symmetry of the problem is such that the m-multipoles (proportional to e irntp ) decouple. 
In the single (barotropic) fluid case, r-mode solutions exist for each individual I = m > 1. 
The main focus in previous work has been on the quadrupole mode (I = 2) since it is 



associated with the fastest growing gravitational- wave instability [32|. In this work we will 
end up studying a large range of scales, including small scale modes with I the order of 100. 
In a neutron star, such perturbations would be relatively local, corresponding to a typical 
lengthscale of 100 m. 

As hinted at in the previous section, we prefer to work with slightly different perturbation 



variables. Experience from other problems involving the two-fluid model 20] suggests that 



it is advantageous to work with the total perturbed momentum flux, defined as 

pU l = pjj[ + p p U l p (22) 
with p = p Q + pp. We take the second variable to be given by the velocity difference; 



.1 jji 



u 



Ul - U{ (23) 



In order to "simplify" the final perturbation equations it is useful to express the frequency 
in the "rotating frame", i.e. work with a defined by 



It is also convenient to introduce 



to + mtt = aVl (24) 

L = l(l + 1) (25) 
Z = 1 - x p - e p (26) 



and 



Z = (27) 

1-Xp 

Finally, we use (as in 20] ) the scaled mutual friction coefficients; 

B' = B'/x p , and B = B/x p (28) 



The perturbation equations can now be obtained by inserting the expected velocity field 
in the equations from the previous section. After some manipulations this leads to the 
following set of the equations to be solved: From p n W^ + p p W p we, first of all, get 



[La - 2m] U% m = —mA(L - 2)[U l - Zu l ]Y t 



(29) 



(Here, and in the following, summation over I > m is implied.) The advantage of working 
with this combination, that represents to total vorticity, is that there are no mutual friction 
terms in the equation. Such terms are associated with the relative flow. To see this we 
consider the combination W p * — W„, leading to 



LoZ 



Xr 



- 2m(l - B') - 2iB{L - m 2 



U% 



= —mA 



[(L - A)x p + 2\Z 



Xr 



2(1 - x p ) | u l Yr + mA^L- ^ U% m 



+ mALB'(Zu l - U^Y™ - 2mB'A(Z + 1 - x^Y™ 
+ iBAL[Z(rd r u l - u l ) - rd r U l + U^Y™ + 2iBA[L -m 2 )(Z + l- x p )u l Y™ 
+ iBLA[rd r U l - 3U l - Z(rd r u l - 3u 1 )} cos 2 BY? 1 

+ iBA [2rd r U l - LU l - Z(2rd r u l - Lu l )\ cos 6 sin OdgY^ (30) 

We will also use the radial components of the Euler equations. From the combination 

[x p rd r 5¥ p + (1 - x^rd^jY^ - 2QU 1 sin^Y,™ = 

= -AQ(1 - x p ) [(x p - Z)rd r u l + 2Zu l - 2U l ] sin ed e Y™ (31) 

while — E r n leads to 



(rd r 5% - rdrS^Y^ - 2fi(l - B')u l sin Qd^ - 2imB£lu L cos 6Y t 

Z 



-Alii f 1 - - > ~* TTl ' 22 T '' 



rd r U l + — U l + (1 - 2x p ) 



x 



1 - 



Xr 



rd r u 



2Z 



-u 



Xr 



-2(l-x p )u l jsm9d e Y l m 

+ B'AQ[rd r U l - Zrd r u l - 2(1 - x p + Z)u l ] smOdeY^ 

+ imBAVt[rd r U l - Zrd r u l + 2(1 - x p + Z)u l ] cos OY^ (32) 
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Finally, it is convenient to work with a "divergence equation" 20] that follows from the 
combination 

sin Ode [sin 9(p n 8 6 n + p p S 6 p )} + d v [p n 8* + p p 8%] — ► 

- L[x p 8¥ p + (1 - x^d^Y^ + 2nU l [Lcos6Y l m + wnOdeXT] = 
= -AQ{1 - x p ) {2U l [L cos 6Y™ + sinfldeY/ 1 ] - Lx p u l [2 cos 6Y™ + sinfldeY™] 

+(L- 2) Zu l sin OdgY™} (33) 

Analogously, we consider the "difference" equation; 

sin Ode [sin 9(8° - £ e n )} + d v [8* - ££) — ► 

L{5V n - ^p)YT + 2fl(l - i3>'[L cos #17" + sin^Y/ 1 ] 

+ 2imQBu l [2 cos BY™ + sin 6d e Y™] = 

= AQL(U l - x p V)[2cos^ m + smOdeY^} 
\L-2)Z 



- AQ 



Xr 



-{U l + u l ) - (L - 2)(1 - x p + 2Z)u' 



sinfl«9 e Y," 



+ Afii3' {-Lf/' + [2(Z + s p ) - (1 -x p )]Ln'} [2cos^ m + sin^YJ" 1 ] 
- AftS' (L — 2)(Z + x p )u l sin #<9 e Y, m 
- 2imttBA(x p + Z)n'[2cos^ m + sin^Y," 1 ] 
- imVtBA[2rd r U l + LC/' + (1 - 2x p - Z)(2rd r u l + Lw z )] cos#Y/ m (34) 

Next, we separate the /-multipoles by means of the standard recurrence relations 

cos6Y l m = Q l+1 Y™ 1 + Q l Y™ 1 (35) 



and 



where 



sin ed e YT = IQ 1+1 Y^ X -{1+ 1)QTO 



(I — m){l + m) 
(2/-l)(2/ + l) 



These relations lead to 



cos^ 6Yr = (Qf +1 + Qf)Yr + Q l+1 Q l+2 Y% 2 + QiQi-^ 



(36) 
(37) 

(38) 



and 



cos 9 sin 6d e Yr = [lQf +1 - (I + l)Qf ] Y™ + lQ l+1 Q l+2 Y% 2 - (/ + l)QiQ;-i^ 2 (39) 
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It should be quite clear at this point that the problem we consider is rather complex, even 
for purely axial modes. One must also be careful, because it is not clear from the outset 
that modes of this particular character exists (as, in principle, rotation would couple the 



axial/toroidal degree of freedom to the polar /spheroidal one 22|). However, in the particular 
case that we are considering the problem simplifies in an almost miraculous fashion. We 
do not expect this level of simplification in a more general situation, leaving the problem 
exceedingly difficult. 

IV. THE UNSTABLE R-MODES 

The r-modes are very special members of the general class of inertial modes, because their 
eigenfunctions "truncate" at I = m (at least in the standard single-fluid setting) making 
their eigenfunctions particularly simple (generic inertial modes involve coupling a number 
of I > m multipoles 22|). Inspired by the fact that this remains true also for co-rotating 
superfluids 20], it makes sense to ask whether it may be the case when a rotational lag is 
present as well. Somewhat to our surprise, it turns out that such simple r-mode solutions 
do, indeed, exist. Moreover, we find that these modes may become unstable once the vortex 
mediated mutual friction is accounted for. 

Assuming that the perturbed velocity fields take the same form as in the co-rotating case, 

i.e. 

, A,r m+1 for / = m 
K={ (40) 
I for I > m 

and noting that Qi =m = 0, the equations from the previous section collapse to two scalar 
relations for the amplitudes. Expressing these in terms of U m and u m , we have 

[(m + 1)<7 - 2 + A(l - x p )(m - l)(m + 2)]U m 

- (1 -x p -£p)Ax p (m- l)(m + 2)w m = (41) 
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and 



- [(m - l)(m + 2) + 2e - m(m + 1)(B' + iB)] AU m 

+ {(1 - e)(m + 1)(7 - 2(1 - B' + iB) + Ai p (m - l)(m + 2) 

{[m(m + 1) - 4]x p + 2}-m(m + l)Ax p (l-£)(i3 / + iB)+2(l-£ n )(B / -ii3)A} 



w m = 



(42) 



where we have used e = e n /x p = e p j (1 — x p ). In the present case, where all the coefficients 
are taken to be constant, it is easy to see that the problem reduces to a quadratic for a. The 
general solution to this quadratic is, however, rather messy and not particularly instructive. 
In order to understand the nature of the solutions, it is better to focus on simplified cases. 

A. Dynamical instability 

Let us first consider the case of vanishing entrainment. In this case we have 



P = 2(m - l)(m + 2) (& - iB) A 2 + 4 [l - (1 + x p ) (B' - iB)] 

+ 2A [(m - l)(m + 2)(x p B' - 1) + (B' - iB) (m 2 + m - 4) + ix p B(m 2 + m + 2)] (44) 



The explicit solutions are, of course, still not transparent. However, if we consider the 
limit of strong mutual friction coupling (TZ —t oo), such that B ~ and B' ~ 1, then we 
obtain the roots 



(m + I) 2 a 2 + P x a + P = 



(43) 



with 



and 





(46) 



with 




(47) 
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These solutions highlight one of the main new results in this paper (and 12] ): For m 3> 1 
we have 

V « (1 + x p ) 2 - 2x p m 2 A , (48) 
showing that we have unstable r-modes (Im a < 0) for 

1 + x D 



> 



where 



y/2x p A 



(49) 



Expressed in terms of the (observed) rotation period P = 2tc/Q of the system, the growth 
timescale for these unstable modes is 

,2 \ -1/2 



mP ( x 



'grow ~ 



27T V 1 + X r 



m | 



(50) 



For m ^> m c (in practice, m > 2m c ), the growth rate is well approximated by 



grow 



2tt \2A/ 



T ~ *? 

' grow ~ 



2 n \l/2 /10- 4 \ 1/2 



0.05 



A 



P 



(51) 



With the same scaling, the critical multipole beyond which the instability is present is 

X n N-l/2 / A N _1/2 



m r rj 300 



0.05 



io- 



(52) 



The presence of this critical value, and the associated emergence of unstable modes, is 
illustrated in Figure [U We see that the instability sets in as two r-modes merge at a critical 
scale represented by m c . This is the characteristic behaviour of a dynamical instability in a 
non-dissipative system [23 ]. 

Accounting for the entrainment obviously complicates the analysis. However, the result 
remains transparent in the strong coupling limit. Again setting B = and B' — 1, we find 
the roots 

a = (w+ 7 1)xp [-(1 + e n ) + (1 - e n )(x p — A) ± (53) 

with 



7 



and 
V 



7 



+ 2(1 - e n ) 



1 E n Bp 



1 - (m 2 + m-3)^ 



-i 



(54) 



7 



A 



+ (l-£n) 



1 - ^(m + 2)(m- 1) 

7 



A 2 (55) 
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FIG. 1: Real (solid lines) and imaginary parts (dashed) of the roots to the dispersion relation (the 
r-modes) in the strong-coupling limit, 1Z = 10 3 , for two illustrative cases, when e p = (thin lines) 
and e p = 0.6 (thick lines). The rotational lag is fixed to A = 5 x 10~ 4 . In each case, the presence 
of a critical value for the azimuthal index m c , beyond which the modes are unstable (the roots 
are complex) is apparent. The results also demonstrate how the entrainment affects the range of 
the instability by shifting m c . Meanwhile, the growth rate of the instability (the magnitude of the 
imaginary part in the unstable regime) is not affected much. 

For m ^> 1 and A <C 1 this is approximately, 

DM (7±xj!A_ where ml = (56) 

7 V m cJ 2x p A7(l-£: n ) 
We see that, in this case there are unstable modes with growth time; 

D / \ / 2 \ -I/ 2 

mr I i n \ / m ' 



^"^Ur^Ur 11 (57) 



For m 3> m c this reduces to 



where we have introduced 



^ P /£p\V^ 1/a 
~ 2tt V2A/ 

1 - e n - £ p _ 
1 - e n m p 



el' 2 (58) 



(59) 

i iip 

Given that a typical range of values for the effective proton mass, m* in a neutron star core 



is expected to be in the range m*/m p « 0.5 — 0.9 28j we conclude that entrainment has a 



minor impact on the growth timescale. An illustration of this result is provided in Figure [H 
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In summary, these resu 
from radio pulsar glitches 



ts show that, for the typical magnitude of rotational lag inferred 
we would have unstable modes with a characteristic horizontal 
length-scale of tens to hundreds of metres. Smaller scale r-modes would be unstable, with 
a growth time as fast as a few rotation periods. As the unstable modes grow extremely 
fast compared to the evolutionary timescale (spin-down, cooling etc) of the system it seems 
reasonable to expect that they may affect any real system that develops the required lag, 
A. It is worth noting that the predicted growth time is much shorter than the current 
observational constraint for the rise of pulsar glitches (tens of seconds [24(), unless the star 
is slowly rotating. Hence, the instability could grow fast enough to serve as trigger for the 
observed events. 



B. Secular instability 

Having established the existence of an instability in the strong-coupling limit, let us 
consider the problem for less "extreme" parameters. It is, of course, straightforward to solve 
(|43|) for given parameter values. A sample of results obtained by considering the problem 
for fixed 7Z, leading to the mode frequencies depending on m, are provided in Figs. [2HH 
These graphs show the behaviour of the r-mode solutions as 7Z decreases from 100 to 2, i.e. 
as we move away from the strong-coupling regime. The results are very interesting. First of 
all, we see that the general trend of a "dynamical" instability (associated with modes with 
a markedly larger imaginary part to their frequency) setting in near the critical value m c 
remains down to 1Z = 10 or so. For smaller values, e.g. 7Z = 2 as in Fig. HJ there is no longer 
a clear change in the imaginary part near m c . We also see that, away from the extreme 
strong coupling limit the dissipative aspect of the mutual friction becomes important. In 
particular, it leads to the dynamical instability no longer being associated with exact mode 
mergers. Rather, the instability sets in at "near misses" in the complex frequency plane. 
This is probably what should be expected. Another aspect of the results was not expected. 
Considering the imaginary parts of the roots in more detail (the right-hand panels of Figs. [2]- 
H|) we see that neither imaginary part changes sign in the displayed interval (as we show 
log I Im a I a sign change would show up as a sharp singularity) . This demonstrates the 
most interesting new result in this paper. As we know that one of the modes is unstable 
beyond m c we must conclude from Figs. [2HH that this mode is, in fact, unstable for all lower 
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values of m, as well. Of course, for smaller m the growth time of the unstable modes is much 
longer. Later we will demonstrate that it is linked to the mutual friction parameters, making 



this a secular instability 
laboratory superffuids 



plausibly related to the Donnelly-Glaberson vortex instability in 



25|-|27|). Particularly interesting may be the fact that this instability 
is not restricted to the small scale r-modes. It is also active for the large scale modes. This 
is interesting since these modes, expecially the m = 2 r-mode, are also secularly unstable 
due to gravitational- wave emission. Basically, our results show that the mutual friction 
may not provide damping of these modes. Rather, it could provide an additional driving 
mechanism for the instability. It may also be, given the strong scaling of the gravitational 
radiation reaction with the modes oscillation frequency (essentially the star's spin rate) that 
the mutual friction driven instability dominates for slowly rotating systems. We will return 
to this question later. 
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FIG. 2: Real (left panel) and imaginary parts (right panel) of the r-modes for 1Z = 100, x p = 0.1 
and a rotational lag A = 5 x 10 -4 . Beyond a critical critical value for the azimuthal index, m c , 
the modes are dynamically unstable. The onset of this instability is associated with (near) merger 
of the real part of the two frequencies. The absence of sign change of the imaginary part of the 
unstable branch shows the presence of a secular instability for smaller values of m. 



Let us first see if we can find approximate solutions that demonstrate the behaviour 
seen in the numerical results. Returning to ( I4ip and ( 142]) we see that the two degrees of 
freedom (represented by U m and u m ) only couple at order A. Since we are interested in 
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FIG. 3: Same as Fig. [2] but for K = 10. 
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FIG. 4: Same as Fig. [2] but for K = 2. 



mode solutions for small A, this suggests that the solutions are either predominantly co- or 
counter-moving, depending on whether U m dominates over u m , or vice versa. Considering 
first the co- moving case, we assume that U m ^ and u m ~ 0, which leads to 

1 



a 



m + 



- [2 - A(l - x p )(m - l)(m + 2)] = a 



(60) 



At this level, the frequency is real and the modes are stable. 
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Meanwhile, in the counter-moving case we have modes with frequency 



a 



_ ) -^{2(1 -B' + iB) - Ax p {m - l)(m + 2) + eA {[m{m + 1) - 4]x p + 2} 

y± Sj \ Tfi ~r 1 J v. 

+ m(m + l)Ax p (l - e){B' + iB) - 2(1 - e n )(B' - iB)A^ (61) 

From this we see that the imaginary part of the frequency is proportional to 

2 + m(m + l)Ax p (l - e) + 2(1 - e n )A 

which suggests that the counter-moving modes tend to be stable, at least in the non- 
entrainment case. When the entrainment is accounted for, these modes may become unsta- 
ble. In the short lengthscale limit, when m 2 x p A ^> 1 (corresponding tom> 10 3 or so for 
typical parameters), an instability is present as long as e n > x p , which is not a particularly 
severe criterion. 

Let us now consider the co-rotating modes at the next level of approximation. We have 
already seen from ( |60l) that, to leading order in A, these modes are marginally stable, with 
frequency approximated by oq. Yet, we know from the numerical results that an instability 
should be present for some range of TZ, certainly 1Z > 2. Hence, we need to estimate the 
solution at the next order in A. To do this, we take 

a = a + a 2 A 2 (62) 

which leads to 

a 2 = ^'^^^I-i] [(m - l)(m + 2) - m(m + 1)(B> + iB)] (63) 

We are (primarily) interested in the imaginary part of the frequency. Since we know that an 
instability exists in the strong-coupling limit we note that this limit corresponds to B' ~ l/x p 
and hence we have (to order B) 

(m — l)(m + 2) sr / w n n ^ / \ 

a2 = 2(m + l) p(1 _ x p)H™ + 1)(2 - x p ) + 2x p ]B (64) 

We learn that these modes are unstable for all values of A. We also see that the growth rate 
scales as 1/B making this a secular instability. The estimate (]64"1) establishes the presence 
of the instability for a wide range of parameters, and provides more detailed insight into the 
dependence on the parameters. 
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FIG. 5: Same as Fig. [2] but for 1Z = 0.5. The inset in the right panel shows the detailed behaviour 
of Im a as function of m. The dashed curve in the right-hand panel corresponds to the approximate 
solution (|63|) . which is an accurate representation for low multipoles. 



It is obviously relevant to establish what the critical value of TZ may be. We can get an 
idea of this by working out where the imaginary part of (|63|) changes sign. This leads to 

TZ\ (m-l)(m + 2) 



(65) 



l + TZ 2 c 2m{m + 1) 

For TZ < TZ C the system should be stable. It is easy to show that this is the case by 
considering the weak-coupling limit. Then we have B' ~ B 2 so the dominant behaviour will 
be 

(m — 1 \ 2 lm -i- 0\ 2 1 

(66) 



(m-l)V+ff 1 
2(m+l) pV p7 i3 



As expected, the modes are always stable in this limit. This is, of course, as expected. 

Returning to the numerical results, the typical behaviour for TZ C < 1Z < 1 is similar to 
that shown in Fig. [5] (which corresponds to 1Z = 0.5). That is, for a given value of TZ the 
instability is present for a range of multipoles, up to a critical value where the modes become 
stable. Above TZ = 1 all modes are unstable. 

We get a complementary view of the new instability by fixing the multipole m and varying 
TZ. This leads to the results shown in Figures [6] and [7J for m = 2 and m = 100, respectively. 
These results demonstrate that the instability exists throughout the TZ > 1 regime, and that 
is extends into the "weak-coupling" regime, as well. This is yet another demonstration of 
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the generic nature of the superfluid r-mode instability. It is notable that the approximation 
fl64|) is excellent for low values of m, see the right-hand panel of Fig. |6l but deteriorates as 
m increases. This is not surprising since, for large values of m one cannot simply neglect 
higher order terms in A as these may be multiplied by factors of m. In essence, the co- 
and counter-moving degrees of freedom are no longer neatly decoupled on shorter (angular) 
scales. Nevertheless, the approximation serves its purpose by illustrating the behaviour in 
an important part of parameter space. 




0.01 0.1 1 10 100 0.01 0.1 1 10 100 

R R 

FIG. 6: Real (left panel) and imaginary parts (right panel) of the r-modes for m = 2, x p = 0.1 and 
a rotational lag A = 5 x 10~ 4 . Unstable modes are present beyond a critical value of 1Z, see inset 
in the right panel. We compare the obtained imaginary part for the modes that become unstable 
to the estimate (|63p (dashed curve in the right-hand panel). In this case, this is clearly a very good 
approximation. The inset in the left panel shows that the real parts of the mode frequencies cross 
at a low value of 1Z, seemingly unrelated to the onset of instability. 



V. ASTROPHYSICAL CONTEXT: ACCOUNTING FOR SHEAR VISCOSITY 

As suggested in [ijj], the superfluid instability discussed in the previous section could 
be relevant for astrophysical neutron stars, in particular glitching pulsars. A "minimum" 
requirement for this to be the case is that the instability grows fast enough to overcome the 
dissipative action of viscosity in neutron star matter. 

A mature neutron star core is sufficiently cold to contain both superfluid neutrons and 
protons (recent evidence suggest that this is the case for a core temperature T < 5 — 9 x 
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FIG. 7: The same as Figure [6j but for m = 100. In this case, (|63p no longer provides a good 
approximation to the imaginary part. 



10 8 K 29_|, |30| . Under these conditions the fluid motion is primarily damped by vortex 
mutual friction (which we have already accounted for, and which drives the instability we 
are considering) and shear viscosity due to electron-electron collisions 31] • Drawing on the 
considerable amount of work that has been done on the gravitational-wave driven r-mode 
instability 32J, we can estimate the shear viscosity damping timescale using an energy- 
integral approach; 

-^modc I rn\ 

r sv = (67) 

where -E mo de is the mode energy and E sv is the shear viscosity damping rate. The damping 
timescale ( IBTj) can be easily calculated using the two-fluid r-mode results of the previous 
sections. However, as we are only interested in a rough estimate we use the result for 
ordinary single- fluid r-modes in a uniform density star [33j] to approximate r sv . This leads 
to 

3 M 1 1.3 x 10 6 /0.05\ 3/2 R 7 6 /2 



A-nr] ee R (2m + 3)(m 



V 



M 



1/2 



(68) 



where Rq = i?/10 6 cm and Mi. 4 = M/1AM & represent the radius and mass of the star, 



respectively. The relevant shear viscosity coefficient (^ e e = Pp^ee) has been taken to be 31| 



1.5 x 10 19 



0.05 



3/2 



irZrr-2 -1 - 



(69) 



where p u = p/10 14 g cm" 3 and T 8 = T/10 8 K. 

We are now in a position where we can compare the growth timescale, r gr0 w, to that 
due to shear-viscosity damping, r sv . Such a comparison is provided in Figure |8] for modes 
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FIG. 8: The dynamical r-mode instability growth and viscous damping timescales r gr0 w and t sv 
(from Eqns. (j57|) and (168p . respectively) as functions of the multipole m for typical pulsar param- 
eters: P = 0.1 s, x p = 0.1 and strong drag 1Z = 10 4 . The stellar mass and radius are fixed at 
the canonical values Mi. 4 = Rq = 1. The shear viscosity damping rate r sv is shown for two core 
temperatures: T = 10 7 K and T = 5 x 10 K. The spin-lag A is fixed to 5 x 10~ 4 and we show 
results for two values of the entrainment parameter, e p = and e p = 0.6. The grey region indicates 
the range of unstable r-modes for T = 10 7 K. 

exhibiting the dynamical instability (in the strong-coupling regime) and typical neutron- 
star parameters. The timescales are shown as functions of the multipole m, and we provide 
results for two choices of the entrainment parameter e p , cf. Fig. [TJ The damping rate r sv 
is shown for two representative core temperatures T = 10 7 K and 5 x 10 7 K, representing 
mature neutron stars. Dynamically unstable superfluid r-modes exist for the range of m 
above the r grow curve but below the r sv curve. The results show that the r grow profile levels 
off (as a function of m) for m > m c and that the instability can overcome the viscous 
damping for a range of scales. We also see that the entrainment has a small effect on the 
asymptotic behaviour of r grow but can significantly affect the critical multipole m c . Finally, 
it is apparent that the range of unstable modes decreases as the star cools. This is an 
interesting observation since glitches are only seen in relatively young pulsars. The results 
in Figure M indicate that shear viscosity would prevent the instability from developing soon 
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after the star has cooled below 10 7 K. It is worth keeping in mind that the core temperature 



may remain above 10 8 K for (at least) the first 10 5 years of a pulsar's life 36]. 

The predicted spin-lag for the dynamical instability to set in also makes a connection 
with pulsar glitches seem plausible. Balancing the mode growth and the viscous damping, 
i.e., setting r grow = r sv , we find the critical spin lag A c above which the instability is active. 
Combining (1S"T]) and (|68|) and setting m = m c we obtain 



m p 



This result was first derived in 



121 ] . As discussed in that paper, the predicted critical lag 
compares well with the available data for large pulsar glitches. This could be an indication 
that large glitches are indeed triggered by the large m r-mode instability. 

The results for the secular instability for smaller values of m are quite similar. In Fig. |9]we 
show results both for m = 2 and m = 100. In each case we see that unstable modes will be 
present above a certain temperature. Keeping in mind that the critical temperature for core 
superfluidity is expected to be below 10 9 K, while a typical glitching pulsar like the Vela is 
expected to have core temperature just above 10 s K, these results make it seem plausible that 
the secular instability can become active in these systems. Whether astrophysical systems 
with 1Z in the required range for the instability exist is not clear at this point. Most work 
has focussed on smaller values, in which case the system would not exhibit the instability 



we have discussed here, but there have been suggestions of larger values [34j . It is also 
worth noting that a stronger mutual friction may help reconcile the discrepancy between 
our theoretical understanding of the r-mode instability with observed astrophysical systems 



35]. 



VI. FINAL REMARKS 



We have provided a more detailed analysis of the superfluid r-modes in a simple uniform 
solid body) differential rotation between the two components. As 



density model with 

already discussed in 12], the results of this analysis are interesting, both conceptually and 
from the point-of-view of astrophysical applications. It is obviously interesting that the 
vortex-mediated mutual friction may lead to the presence of an instability. However, this is 
perhaps not surprising given the wealth of results relating to instabilities triggering superfluid 
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FIG. 9: The secular r-mode instability growth and viscous damping timescales r grow and r sv (ob- 
tained from the numerical solution to (|43p and (|68p . respectively) as functions of 7Z for two mul- 
tipoles: m = 2 (left panel) and m = 100 (right panel). The other parameters are the same as in 
Fig. S P = 0.1 s, rE p = 0.1 and M 1A = R 6 = 1. The spin-lag A is fixed to 5 x 10~ 4 . The shear 
viscosity damping rate r sv is shown for various temperatures as indicated in each panel (dashed 
horizontal lines). The r-mode growth time is short enough for an instability to be present for a 
range of 1Z for temperatures characteristic of young neutron stars. 



turbulence 
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381 ] and previous results in the context of pulsar precession |39|, |40( , but the 
present analysis provides the first demonstration of this kind of instability for global mode 
oscillations. We are also breaking new ground by considering perturbations for background 
configurations with differential rotation. It is notable that, as soon as we allow for this extra 
degree of freedom the problem becomes richer. 

A particularly interesting aspect of the results is the presence of both secular and dynam- 
ical instability behaviour in (more or less) distinct parts of parameter space. The behaviour 
that we have unveiled is summarized in the phase plane in Figure [TUJ that shows the unstable 
region in the m — 1Z plane. As far as we are aware, the present analysis represents the first 
detailed study of a problem that has secularly unstable modes entering a regime where they 
become dynamically unstable. The associated behaviour is, in many way, predictable. For 
example, instead of having modes becoming dynamically unstable as real-valued frequency 
pairs merge and form complex conjugate pairs, we now have dynamical instability behaviour 
associated with "near misses" in the complex frequency plane. Our results suggest that dy- 
namical instability only results provided that the damping of the modes is not too large (in 



24 



our case 1Z > 10 or so). To improve our understanding further, we need to consider the 
necessary and sufficient criteria for the superfluid instability to operate. This is a challenging 
problem but it seems likely that one could make progress by adding the perturbed mutual 



friction force to the results in 



43|. 



An interesting analogous problem concerns the I = m = 2 bar-mode in rotating neutron 
stars. This mode is known to become secularly unstable due to the emission of gravitational 
waves well before (along a sequence with increasing degree of differential rotation) it becomes 
dynamically unstable, see 23) for a review and 4jj for recent work. The behaviour of that 
instability may be quite similar to that of the superfluid problem we have discussed here. 
If this is the case, then a key question concerns whether there always is a dramatic change 
in gravitational- wave emission rate before and after the critical parameter for onset of the 
dynamical instability is reached. In the bar-mode problem there is also another class of 
instability, commonly referred to as the low T/W instability j^J. We have no evidence for 
an analogous instability in the present problem. 

Even though we have not discussed possible astrophysical repercussions in any detail, it 
is clear that the potential link with the unresolved problem of radio pulsar glitches provides 
strong motivation for further work on this problem. Of course, we cannot at this point 
really tell how relevant the new r-mode instability is in this context. In order to act in 
an astrophysical system, the instability must be robust enough to remain once we account 
for the star's magnetic field and the elastic crust. These features are likely to affect the 
instability considerably, but more work is required if we want to quantify what the effects 
may be. The two-stream instability mechanism is sufficiently generic that it would be 
remarkable if it would cease to operate in more complex settings, but it could be that the 
instability threshold moves out of reach for a real system. At this point, we cannot say. 
There is, however, fresh evidence supporting the existence of short wavelen gth superfluid 
instabilities in the presence of a magnetic field and/or an elastic crust, see 44], |45|. The 
exact relation between these instabilities and our r-mode instability is still unclear, apart 
from the fact that both require vortex mutual friction to operate. The secular instability 
behaviour is obviously also worthy of further attention. In particular, since it is relevant 
also for large scale modes, including the lowest multipoles that are the most important from 
the gravitational- wave point-of-view. The problem is clearly extremely interesting and well 
worth returning to in the future. 
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changes from secular to dynamical instability on the strong coupling regime (upper panel), and 
the estimated critical drag 7Z C from (165j) where the secular instability sets in for low m-multipoles 
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